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ABSTRACT 

One well-known way to constrain the hydrogen neutral fraction, xhi, of the high-redshift in- 
tergalactic medium (IGM) is through the shape of the red damping wing of the Lya absorption 
line. We examine this method's effectiveness in light of recent models showing that the IGM 
neutral fraction is highly inhomogeneous on large scales during reionization. Using both an- 
alytic models and "semi-numeric" simulations, we show that the "picket-fence" absorption 
typical in reionization models introduces both scatter and a systematic bias to the measure- 
ment of xhi- In particular, we show that simple fits to the damping wing tend to overestimate 
the true neutral fraction in a partially ionized universe, with a fractional error of ~ 30% near 
the middle of reionization. This bias is generic to any inhomogeneous model. However, the 
bias is reduced and can even underestimate xm if the observational sample only probes a 
subset of the entire halo population, such as quasars with large HII regions. We also find that 
the damping wing absorption profile is generally steeper than one would naively expect in 
a homogeneously ionized universe. The profile steepens and the sightline-to-sightline scatter 
increases as reionization progresses. Of course, the bias and scatter also depend on xhi and 
so can, at least in principle, be used to constrain it. Damping wing constraints must therefore 
be interpreted by comparison to theoretical models of inhomogeneous reionization. 
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1 INTRODUCTION 

The reionization of hydrogen in the intergalactic medium (IGM) 
is a landmark event in the early history of structure formation, be- 
cause it defines the moment at which galaxies (and black holes) 
affected every baryon in the Universe. As such, it has received a 
great deal of attention - both observationally and theoretically - in 
the past several years. U nfortunately, the existing observational ev- 
idence is enigmatic (see lFan et al Hood for a recent review). Elec- 
tron scattering of cosmic microwave background photons implies 
that re ionization occurr ed at z ~ 10, albeit with a large uncer- 
tainty JPage et alj[2007h . On the other hand, Lya forest spectra of 
quasars at z ~ 6 show some evidence for a r apid transition i n 
the globally-averaged neutral fraction, Shi (e.g., Fan et al. 2006). 
Howev er the Lya absorption is so saturated in the lGunn & Peterson! 
dl965h trough (with optical depth tgp > 10 xhi) t h at constraints 



derive d from that spectral region l Fan et alj [20061; iMaselli et al 
2007 ) are difficult to interpret (e.g. Ilidz et allbOOfjUBecker et al 
20071) . 



Another probe is the red damping wing of the IGM Lya ab- 
sorption: the line is so saturated at these redshifts that even photons 
that are emitted redward of the Lya resonance can suffer significant 
absorption from the strong damping wings of that transition. This 
has a number of consequences for high-redshift observations. 

For example, surveys that search for high-z galaxies through 
their Lya emission lines will find fewer and fewe r galax- 
ies as the IGM becomes more and more neutral jHaimanl 



2002; ISantosI 12004). although galaxy clustering strongly mod- 
erates this decline ^Furlanetto et al.l 2004 120061: IMcOuinn et al.1 
1 20071 : iMesinger & Furlanetto! 20071 ). Such su rveys have now 
detec te d object s at z ~ 6. 5 -9 ( e.g., iKashikawa et al.l 
120061 : live et all 120061 ; IStark et alj l2007t). but their implica- 
tions for reionization are unclear jMalhotra & Rhoadsl |2004 



lHaiman & Cen 2005; Malhotra & Rhoads 2006; Kashikawa et al 



2006; Dawson et al. 2007; Diikstra et alj 120071 : IMcOuinn et al 
120071 : [Mesinger & Furlanettdl 2007b!) . 

The evolution of galaxy abundances and clustering measures 
the damping wing absorption in a statistical sense, but even more 
information can potentially be gleaned from the damping wi ng ab- 
sorption profiles in individual objects ( Miralda-Escude 1998). For 
the galaxies described above, this information is difficult to extract 
because of their faintness and the com plicated origins of their Lya 
emission lines ( IMcOuinn et alj2007l) . 

However, high signal-to-noise spectra of bright objects could 
be extremely helpful. If the damping wing profile from IGM ab- 
sorption can be isolated from these spectra, this would provide de- 
tailed information on the neutral gas along each particular line of 
sight (LOS) - rather than the statistical information available from 
most other probes. This is very useful, as reionization is expected 
to be highly inhomogeneous. 

There are two candidates for such high signal-to-noise spectra 
at high-redshifts: quasars and gamma- ray bursts (GRBs). Quasars 
present several cha llenges: complicated intrinsic spectra, b iased 
IGM environments (Barka na & Loebl l2004; Lidz et al. 2007]), and 
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large HII regions (which significantly weaken the damping wing 
absorption redward of the quasar Lya line, and ca n necessitate de- 
tailed spectral analysis of t he blue side of the line: lMadau & Reesl 
2000; Mesi nger et al]2004h . Nevertheless, there have already been 
two claims of damping wing detections in high-redshift spectra, 
both using quasars from th e Sloan Digital Sky Survey (SDSS). 
iMesing er & Haiman (2004) detected a xm > 0.2 damping wing 
through the decreased fluctuations in the total Lya optical depth 
near the edge of the HII region surrounding J1030+0524 {zs = 
6.28). Similarly, by simulating the optical depth distributions blue- 
ward of th e Lya line center and comp aring them with deep ob- 
servations, Mesinger & Haiman (2007) detected the presence of 
a xm > 0.033 damping wing in the spectra of J1030+0524 
and J1623+3112 (zs = 6.22). The maximum likelihood was at 
xm = 1 for both quasars. 

The second set of candidates, GRBs, have fewer obstacles 
to overcome. Long-duration GRBs are believed to be remnants 
of massive stars (and so trace the bulk of the star formation, 
which probably occurs in lower-mass halos with more "typical" 
IGM environments), and their afterglows have extrem ely sim- 
ple power-law intrinsic spectra (see, e.g., iPiranl 120051 for a re- 
view). The event rates at high redshifts may be quite high, and 
cosmological time-d ilation helps to ident i fy the sources when 
they are still bright jBro mm & Loeb ll2002l; ICiardi & Loebll200Cl : 
lLamb & ReicharfeoOOt IMesinger et "ai] |2005h . As a result, there 
is a great deal of optimism in the literatur e regarding their poten- 
tial for damping-wing measurements (e.g. Miralda-Escudel ["l998l : 
iBarkana & Loebll2004h . The highest-redshift GRB afterglow ob- 
served so far (at z ~ 6.3), has already been used to con- 
strain the global neutral fraction at that time dKawai et ai] |2006l ; 
iTotani et al.ll2006h . Unfortunately, this object illustrates the major 
difficulty with the red dampi ng wing test for GRBs : intrinsic ab- 
sorption in the host galaxies dMiralda-Escudel fl998h . Most GRBs 
are now known to have large columns of asso ciated neutral hy- 
drogen dVreeswiik etafl .2004 ; Chen et al ] |2005l). Roughly 20% of 
well-studied objects have 7V H i < 10 20 cm" 2 dChen et al.l2007l) . al- 
though nearly all of the objects in this sample are at z < 6. 

The z ~ 6.3 GRB d oes appear to hav e intrinsic absorption 
with/V H i ~ 10 21 - 6 cm" 2 dTotani et alj2006h . which makes it diffi- 
cult to constrain the IGM absorption. In principle, it is still possible 
because isolated HI absorbers have different spectral profiles than 
the IGM (with the optical depth inversely proportional to the wave- 
length offset squared for isolated absorbers, and to the wavelength 
offset itself for the IGM). The two sour ces can then be sepa rated by 
looking at the shape of the absorption. ITotani et alj d2006t) found a 
best fit with xm = and estimated that xm < 0.17 (0.60) at 68% 
(95%) confidence. Better constraints will require faster followup 
(when the afterglow is brighter) and systems with less intrinsic ab- 
sorption. 

To date, the red damping wing test has generally been as- 
sumed to be simple and straightforward. It is usually argued that 
the absorption is sensitive to a large path length in the IGM, so 
that small-scale dumpiness can be ignored and that th e ionized 
fracti on can be taken to be uniform (for an exception, seejgarkana 
2002). However, most models of reionization have much more in- 
homogeneous distributions of neutral and ionized gas, with dis- 
crete HII regions surrounding clusters o f galaxies, and a sea of 
nearly neutral gas sepa rating them (e.g., lArons&Wingertlll972l ; 
Shapir o & Girouxl 19871) . Such a picture is inevitable when hot stars 
ionize the gas. Moreover, the most recent models show that the ion- 
ized bubbles can become quite large even relativel y early in reion- 
ization, with sizes > 10 Mpc when Xm ~ 0.5 dFurlanetto et al.l 



2004 1 120061; llliev et alj 120061: IZahn et alj |2007| ; iMcOuinn et ail 
20071 : IMesinger & Furlanettoll2007al) . 

Because the damping wing is sensitive to fluctuations on Mpc 
scales, it is actually not a good approximation to take the IGM ion- 
ized fraction to be constant. In this paper, we will examine whether 
(and how) the damping wing can actually be used to constrain the 
reionization process. We summarize the basic physics of the line in 
fj2] We then examine a series of toy models of the "picket-fence" 
absorption typical of the IGM during reionization in ij3] In partic- 
ular, we show that interpreting measurements with the naive view 
of a uniform IGM is not only subject to significant scatter (from 
the different networks of ionized bubbles intersected along differ- 
ent lines of sight) but also a substantial systematic bias. In Sj4] we 
describe the "semi-numeric" simulations used to generate our main 
results, which we present in §5\ This more detailed picture confirms 
that scatter between different lines of sight and bias relative to the 
naive view will be critical in interpreting any observed sources. Fi- 
nally, we conclude in Sj6] 

When this project was nearing co mpletion, we learned of a 
similar effort by IMcOuinn et alj d2007h and refer the reader there 
for a complementary discussion. 

In our numerical calculations, we assume a cosmology with 
Q m = 0.26, Ov = 0.74, Q b = 0.044, H = 100/t km s" 1 Mpc" 1 
(with h = 0.74), n = 0. 95, and cr 8 = 0.8, consistent with the 
most recent measurements dSpergel et ai]|2007l) . Unless otherwise 
specified, we use comoving units for all distances. 



2 THE Lya DAMPING WING 

We compute the total Lya optical depth at an observed wavelength 
Aobs = A a (l + z) along a line of sight (LOS) centered on a halo 
at zs. We do this by summing the damping wing optical depth, td, 
contribution from each neutral hydrogen patch (extending from zm 
to z e i for the ith patch, w ith Zbi > z e Q encounte red along the LOS, 
using the approximation dMiralda-Escudlll998h : 
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where top « 7.16 x 10 5 [(l + z s )/10] 3/2 is the lGunn & Petersonl 
( 1965) optical depth of the IGM in our cosmology, R a — 
A/(47n/ a ), A = 6.25 x 10 s s" 1 is the decay constant for the Lya 
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This expression is only valid far from line center, but that is accept- 
able because the optical depth is so large (and therefore unmeasur- 
able) at line center anyway. It also assumes fi m (z) = 1, which is 
an excellent approximation at the high redshifts relevant here. 

For the remainder of this paper, we will assume that xm(i) = 
1: in other words, the neutral patches between ionized zones are 
indeed fully neutral. This is an excellent approximation in numer- 
ical simulations of reionization by hot stars (although less good if 
X-rays contribute substantially). 

For the analytic calculations in the following section, it is use- 
ful to approximate I(x) by its asymptotic limit when \x — 1| -C 1. 
In that limit, equation l[T) can be written 
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t d (z) 
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where Rbi and R e i are the comoving distances from redshift z to 
redshifts zu and z e i, the beginning and end of the ith neutral patch 
in redshift space. Note that we have further assumed z e i ~ zu ~ z 
(which is an excellent approximation in most of the cases of inter- 
est). Although this asymptotic form is not accurate enough for de- 
tailed calculations or inferences from observations, it contains all of 
the essential features of the damping wing solution and so is useful 
to understand the source of many of the effects we will describe. 
We will also use wavelength units, AA b s = A b s — Aa(l + Zu)- 

The R^ 1 ~ (A ba/AA bs) decline of equation ([3} is much 
gentler than a damped Lya absorber (DLA) at the same location 
(which falls off like [A bs/AA b s ] 2 ). This is because of the large 
sizes of the IGM damping wing absorbers: at large wavelength off- 
sets, a l onger path length is a ble to contribute, which moderates the 
decline dMiralda-Escudell 19981) . It is this property that (one hopes) 
will allow us to distinguish IGM absorption from neutral gas within 
the host galaxies of GRBs, for example. 

A common alternative approach to ours (where we have ex- 
plicitly broken up the LOS into ionized and neutral patches) is to 
assume that the damping wing averages over a sufficiently long 
path length so that the sum over the patches can be replaced by an 
average neutral fraction, xd' 
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or in the asymptotic limit of I(x), in analogy with equation {3), 
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where z b i and Rbi denote the edge of the closest neutral gas and 
z e and R e denote the largest distance at which neutral gas sits (the 
result is, however, quite insensitive to z e , as long as z e ^ z b \). In 
this simple picture, the spectral profile of the absorption is well- 
known. 

The basic measurement is then to extract Zbi (or Rbi) and 
xd from a fit to td(z). It is commonly assumed in the literature 
that xd will be an accurate and unbiased estimator of xm\ in the 
remainder of this paper, we will critically examine these expecta- 
tions. Of course, in principle a much more sophisticated fit may be 
performed with many Rm and R e i. However, in practice the de- 
pendence on any individual element (other than Rbi) is small, and 
adding more parameters will rapidly weaken constraints from the 
fit. 



3 SOME ILLUSTRATIVE TOY MODELS 

We begin by constructing a series of toy models, using the approx- 
imate form of the damping wing optical depth in equation {3j, that 
will show the basic features of the measurement. For this simple 
case, we will define the apparent neutral fraction xd by equat- 
ing the right hand sides of equations l[3} and l[5} and taking the 
R e — > oo limit in the latter: 



X D = Rbi ( ^JCR&i 1 - Rei 1 ) 



(6) 



where we have assumed that the only second parameter measurable 
from the absorption profile is the location of the near edge of neu- 



tral gas, i?i,iQwe will also assume that the summation extends to 
infinity, although this does not affect our conclusions. 

3.1 Bias in the Inferred Neutral Fraction 

To begin, we consider a model where the IGM is divided into ion- 
ized and neutral patches of fixed lengths Rb and fRb, respectively, 
where / = xhi/(1 — Shi) ensures that the mean neutral fraction 
along the given line of sight is Shi- For all of our toy models, we 
will assume that the source sits in the middle of its host bubble, so 
Rbi = Rb /2; keep in mind, however, that this does not affect our 
results because we have assumed that Rbi is independently mea- 
sured from the data. In reality, more sophisticated techniques ma y 
be needed to constrain Rbi (e.g. Mesinger & Haiman 2004. [200% . 
Equation ((6} then becomes 
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The short-dashed line in Figure [TJ shows the difference (xd — 
Xm) for this model as a function of the true average neutral frac- 
tion, xm- It is obviously not a particularly good estimator. The error 
peaks at ~ 0.3 when xm ~ 0.5 (although note that the fractional 
error continues to increase as Shi — > 0). We see that xd overest- 
imates the neutral fraction, because the nearest material dominates 
the absorption. In the "true" model, this material is always fully 
neutral and so absorbs considerably more radiation than would be 
expected in a simple, uniformly ionized model. 

The thin dotted curves in Figure [TJ show explicitly that this 
bias does not depend on our use of the approximate equation (|3jl. 
Here we take the full expression for td{z) (i.e. we estimate xd 
by equating the right hand sides of equations[TJand|4j but consider 
the same sequence of ionized and neutral regions. The three curves 
assume Rb = 1, 10, and 30 Mpc, from bottom to top. All also 
only include neutral patches at z > 6 (with the source at z s — 9). 
Obviously, the 1/R model is an excellent estimate. Interestingly, 
the bias is only a weak function of the actual size of the patches, 
although it does increase slowly with Rb- This suggests that the 
bias can probably be studied fairly robustly. 

The solid line represents a slightly more sophisticated model. 
Here we keep the ionized patches at a fixed size Rb but draw the 
length of each neutral patch from a uniform distribution over the 
range [0, 2 fRb]; thus the mean neutral fraction is still xm, but 
different lines of sight construct it in different ways. In this case, 
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The solid curve shows the bias in the neutral fraction measurement 
for such a model. It is typically about half that of the model with 
fixed path lengths, although the two converge at small xm- The 
bias is smaller in this case because the shorter neutral path lengths 
decrease the apparent absorption by a larger factor than longer seg- 
ments increase it, and the absorption is weighted more heavily to 
nearby gas. 



1 Note that independently measuring Rbi from the spectrum becomes non- 
trivial for large Rbi, as the resonant absorption from resid ual HI within 
Rbi can become strong enough to wipe out detectable flux I Mesinger et al. 
l2004tlBolton~& Haehnelt 2007). 
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Figure 1. The damping wing bias that results from assuming a constant neu- 
tral fraction throughout the IGM. Each curve shows the difference between 
the inferred and true neutral fractions. The short-dashed curve assumes that 
ionized and neutral patches have fixed sizes. The thin dotted curves as- 
sume the same but take the full damping wing profile, for Rb = 1, 10, 
and 30 Mpc, from bottom to top. The solid and long-dashed curves fix Rb 
but assume respectively that the neutral patches are uniformly distributed in 
size between [0, 2 fRb] and exponentially distributed (with the same expec- 
tation value). The dot-dashed curve assumes that both ionized and neutral 
patches are exponentially distributed. 

Next, we take a model where the ionized patches have a fixed 
size Rt but draw the length of each neutral patch from an expo- 
nential distribution, with expectation value fRb- In this case, the 
probability distribution of u = R e k / Rb is 

f- k exp {-/->-(& -1/2)]} 
PekW " (fc-1)] u-Ofc-1/2)]- 1 

u>(fc-l/2), (10) 

and zero elsewhere. The minimum of u is set by the (k — 1) ion- 
ized bubbles that precede this edge, plus the Ri/2 region imme- 
diately surrounding the source. Because the ionized patches have 
fixed size, Pb(h+i) («) =Pek{u~ 1). 

The long-dashed line in Figure [TJ shows the resulting bias. 
Again, it is much reduced from the case with fixed neutral patch 
sizes and is similar in magnitude to the uniformly distributed case. 
However, the bias in this model tends to be larger at small xm and 
smaller at large xm- This is because the long tail of thick absorbers 
allowed in the exponential model is quite significant when most 
absorbers are narrow. However, when most absorbers are thick (at 
high xui), the larger probability of narrow absorbers tends to de- 
crease the bias. 

Our final toy model allows both the sizes of the neutral patches 
and the ionized patches to be exponentially distributed, with expec- 
tation values fRb and Rb, respectively. In this case, the distribu- 
tions p e fc (u) can also be written analytically, but there is no simple 
pattern with k as in equation JlOb . We therefore simply present the 
resulting bias as the dot-dashed curve in Figure [TJ It is nearly the 
same as the model with fixed ionized patch size at xm ;^0.6: in 
this regime, the neutral patches are much thicker and so their scat- 



Figure 2. Fractional scatter in neutral fraction measurements using the 
damping wing, including only the contribution from (he first neutral region. 
The solid and long-dashed curves assume that the width of the region is 
uniformly and exponentially distributed, respectively (as in Fig.fT}. 

ter dominates. On the other hand, at small neutral fractions, the bias 
is much smaller in this model, because the ionized bubbles become 
on average larger than their neutral neighbors 

The different biases between our toy models illustrated in Fig- 
ure [TJ show that the bias does carry some information on the dis- 
tribution of neutral and ionized gas. In principle, an independent 
measurement of xta would then allow additional constraints on the 
reionization morphology. However, this is likely to be a difficult 
game, because the differences are modest in the realistic models. 



3.2 Scatter in the Measurements 

Our last three models draw path lengths from different distribu- 
tions; in addition to the bias, they will also therefore have scatter in- 
trinsic to the measurement of xo- This is more difficult to estimate 
analytically, because there is significant covariance between the lo- 
cations of the different neutral patches (even without covariance in 
their individual lengths, the ith neutral patch must occur before the 
[i + l]th). For a simple estimate, however, we note that most of 
the absorption (on average ~ 80% in our toy models) comes from 
the first region, so we take the variance of the first term in the sum 
in equation l|6}. Again, more sophisticated treatments are possible 
but not necessary given our access to simulations that include many 
more effects than we can hope to add to analytic models. 

Figure|2]shows the standard deviation in these measurements, 
normalized to the true ionized fraction, for two of our models from 
SjXTJ The solid line assumes that the sizes of the ionized patches 
are fixed but that the sizes of the uniform patches are uniformly dis- 
tributed; the short-dashed line assumes that the latter are exponen- 
tially distributed. Here we see that the fractional scatter increases as 



2 Note that the bias in this model appears to become negative at 
SHI ^5 0-04; however, this is a numerical artifact of our approximations. 
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Figure 3. Damping wing absorption profiles, as a function of fractional 
wavelength offset from the source (at redshift z§). The thick curves show 
the "true" absorption profiles for Shi = 0.9, 0.5, and 0.1 (for the dashed, 
solid, and dotted curves, respectively). Note that the two dashed curves 
overlap and are practically indistinguishable. The corresponding thin curves 
show the absorption profiles for uniformly ionized IGM normalized to the 
same transmission at zg. The dot-dashed curve shows the profile of a DLA, 
normalized to the same transmission as the xjji = 0.1 curves at zg. 



Shi — » (although, just as with the bias, the actual value of a xml 
peaks at Shi ~ 0.5). 

Interestingly, the scatter is larger for the exponentially dis- 
tributed model - unlike the bias. This is because of the long tail 
allowed by the exponential distribution which becomes especially 
important at small cehi: the variance of our uniform distribution is 
f 2 R%/3, while the variance of the exponential distribution is f 2 R 2 . 
Obviously, interpreting any observations in detail will require care- 
ful modeling of the underlying distribution. 

As with the bias, the scatter evolves throughout reionization. 
However, it can be measured even without an independent estimate 
of xhi and so can itself be used for further constraints: a large dis- 
persion in the measured x d is indicative of the final stages of reion- 
ization, for example. 



3.3 The Absorption Profile 

The final question we can address with our toy model is how the 
"picket fence" absorption characteristic of inhomogeneous reion- 
ization affects the damping wing absorption profile as a function 
of wavelength. Of course, in our more realistic models that allow 
a range of absorber sizes there will be a corresponding range of 
profiles, and with the large number of absorbers that are allowed 
it is not obvious how one would fit the results or even parameter- 
ize the possibilities. We therefore focus on the simplest model, in 
which the ionized and neutral regions have fixed sizes Rt and fRb, 
respectively. 

Figure [3] shows several example profiles as a function of the 
fractional wavelength offset from the source Lya line center (at a 
redshift zs)- The thick curves are taken from our toy model; the 
dashed, solid, and dotted curves take xm = 0.9, 0.5, and 0.1, 



Figure 4. Residuals between "true" absorption profiles and those from a 
uniformly ionized medium with X£> , as a function of wavelength offset from 
the source. In (a), we set R), = 10 Mpc and vary the IGM ionized fractions. 
In (b), we fix Sjjl = 0.5 and vary the bubble sizes. 



respectively. The thin curves show the corresponding profiles for 
a uniformly ionized IGM, beginning the same distance from the 
galaxy, and with an assumed neutral fraction xd- Thus, they are 
normalized to have the same transmission as the "true" curves at 
the redshift of the galaxy. 

The profiles are nearly identical when Shi is large, but in the 
middle and end stages of reionization the toy model has steeper 
absorption than a uniformly ionized IGM, allowing slightly more 
transmission redward of the source wavelength. This is not surpris- 
ing: as described above, the gentle decline of the damping wing 
occurs because longer columns contribute to the absorption of pho- 
tons far to the red. In the picket-fence picture, photons far from zs 
are sensitive to such a large column that they average over many 
ionized patches. This also explains why the effect becomes more 
and more important at smaller xm, as more and more of the rele- 
vant absorbing gas is absent when the neutral patches are narrow 
and widely separated. 

One possible worry is confusion of the damping-wing absorp- 
tion with DLAs in the host galaxies of GRBs; separating the two 
sources of ab sorption requires th at they have significantly different 
profiles (e.g.. lTotani et alfeoOot) . The dot-dashed curve in Figure|3] 
shows a DLA profile with t(z s ) normalized to the optical depth 
in our picket-fence model with xm = 0.1, where the true pro- 
file is steepest. In this case, the picket-fence model is about mid- 
way between the DLA and IGM expectations, so we would ex- 
pect that clearly separating IGM and DLA absorption will become 
significantly more difficult toward the end of reionization. How- 
ever, when Shi > 0.25, the picket-fence absorption is much closer 
to the naively expected IGM behavior than to the DLA behavior. 
(Of course, if the DLA is centered at zs, it will obscure much more 
of the line profile - but this shows that they can be distinguished, at 
least in principle.) 

Figure [4] shows the differences in more detail. We plot T p f — 
T u , where T p { is the transmission in the picket-fence model and 
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T u is the transmission for a uniformly ionized IGM, normalized to 
the same optical depth at zs- Panel (a) shows the residuals for sev- 
eral different ionized fractions with Rb held constant (as in Fig.O, 
while panel (b) varies Rb while holding Shi constant. 

The deviations are typically at most a few percent, with the 
biggest differences when (z — zs) /(1 + zs) < 0.02. Of course, this 
is also the region most likely to be contaminated by an absorber in 
the host galaxy, so it is not clear how well this part of the deviation 
can be detected. At redder wavelengths, the differences are < 2%, 
so they will require extremely high signal-to-noise spectra to detect 
them cleanly. 

The differences at small offsets from zs are much larger for 
smaller bubbles, even though the bias in the estimated neutral frac- 
tion is nearly independent of Rb (see Fig.[]}. This is because, when 
the first neutral patch is large, most of these photons are efficiently 
absorbed inside of it. When the patch is small, the effective column 
decreases relatively rapidly. Note that the strongest differences in 
the profiles, especially at small wavelength offsets, are due to Rb- 
Thus most of the variations in the spectral shape will help to pin 
down the bubble size (which we have assumed can be measured 
independently). This suggests that it will be difficult to use varia- 
tions in the shape to measure xd more accurately, at least not in 
any straightforward fashion. 

It is obvious from this section that the damping wing spec- 
tra contain more information than just the location of the nearest 
neutral gas, Rbi, and xd (which we have assumed to be the two 
observables). However, it is not clear whether higher-order differ- 
ences can be measured in practice, because of the finite signal-to- 
noise to be expected from these sources and because of interven- 
ing absorption in the host galaxy itself. This is especially true be- 
cause the number of extra parameters required is large: for exam- 
ple, we have found that including only the first neutral region leads 
to residuals of similar magnitude to those in Figure [4] (though with 
the opposite sign, because such a procedure always underestimates 
the total amount of absorption). Thus an accurate fit would require 
adding the contributions from many neutral patches, each with un- 
known location and width. 

Moreover, there will of course be scatter in the profiles at a 
given (xd, Rbi) because of scatter in the sizes of ionized and neu- 
tral patches. In the following we will therefore take the simple- 
minded viewpoint that only these two quantities can be measured, 
although we will use our simulations to measure the dispersion in 
the profiles. We defer a more detailed investigation of parameter 
estimation to future work. 



of the linear density field using excursion-set theory, with mass 
scales spaced as AM/M = 1.2. Note that we are able to re- 
solve halos with masses less than a factor of two from the cool- 
ing mass likely to be pertinent mid-reionizatio n, corresponding 
to gas with a temperature of T ~ 10 4 K (e.g. Ef stathioul 1 1 9921 ; 
iThoul & Weinberg 199dJGnedinll2000l : IShapiro et al.lll994h . Halo 
locations are then adjusted using first-order Lagrangian perturba- 
tion theory. The resulting halo field matches both the mass function 
and statistical clustering pro perties of halos in N-body simulations 
dMesinger & Furlanettoll2007al) . 

In constructing the ionization field, the IGM is modeled as 
a two-phase medium, comprised of fully ionized and fully neu- 
tral regions (this is a fairly accurate assumption in the context of 
damping wing absorption before the end of reionization, unless the 
X-ray background is rather strong). Using the same halo field at 
z = 9, we generate ionization fields corresponding to different val- 
ues of Shi by varying a single effi ciency parameter, C> again us- 
ing the excursion-set approach (c.f. [Mesinger & Furlanettoi l2007a; 
iFurlanetto et alj|2004l) . 

This semi-numeric method is thus ideally suited to the damp- 
ing wing problem, because we are able to "resolve" relatively small 
halos and simultaneously sample a large, representative volume of 
ionized bubbles. 

Our procedure f o r comp uting td is fully described in 
iMesinger & Furlanetto! d2007bl) . with a couple of minor differ- 
ences noted below; no te that similar results were also obtained by 
iMcOuinn et al.l {2007) using full radiative transfer simulations. We 
use equation {T} to calculate the Lya optical depth for each neu- 
tral hydrogen patch, summing the contributions of patches along 
the LOS for 200 Mpc or until the first neutral patch is encoun- 
tered, whichever comes lastQ wrapping around the simulation box 
if needed. We construct distributions of td for halo mass scales in 
the range 2.5 x 10 9 - 2.4 x 10 10 M Q and for various ionization 
topologies (i.e. Shi). We make sure to process LOSs from every 
halo of a particular mass scale, cycling through the halo list un- 
til each mass scale underg oes a minimum of 3 x 10 4 such Monte 
Carlo realizations. Unlike Mesinger & Furlanetto] d2007bl) . we do 
not include the peculiar velocities of halos in estimating tdQ 



5 RESULTS 

In this section, we use the semi-numeric simulations to revisit the 
issues of bias and scatter raised in Sj3] We begin by illustrating the 



4 SEMI-NUMERICAL SIMULATIONS 

In the remainder of the paper, we will use more reliable and 
physically-motivated calculations of inhomogeneous reionization 
that incorporate the full geometry of the IGM to examine the same 
issues of bias and scatter in damping wing measurements. We use 
an excursion-set approach combined with first-order Lagrangian 
perturbation theory to efficiently generate density, velocity, halo, 
and ionizati on fields at z = 9. This "semi-n umerical" simulation is 
presented in lMesinger & Furlanetto! {2007a), to which we refer the 
reader for details. A similar halo -finding scheme has also been pre- 
sented bv lBond & Myers! {l996) and a similar scheme to generate 
ionization fields has been presented bv lZahn et alj d2007l) . 

Our simulation box is 250 Mpc on a side, with the final den- 
sity, velocity and ionization fields having grid cell sizes of 0.5 Mpc. 
Halos with a total mass M > 2.2 x 10 s M are filtered out 



3 This number was chosen experimentally in order to ensure convergence 
of the m distributions at the mass scales and neutral fractions studied in 
this work. 

4 Ignoring velocities simplifies the analysis, since we are guaranteed not 
to have halos whose Lyce line centers have Doppler shifted into resonance 
in the neutral IGM. For such objects, the damping wing could still be re- 
covered by looking further redward in their spectra. The same fundamental 
quantities (especially an analog of xr>) could still be measured from such 
sources, but they would need to be re-parameterized. For the purposes of 
our statistical analysis it is useful to compare absorption at a single redshift 
offset across all objects, which we take to be the line center, 2g . We have 
verified that including halo peculiar velocities mainly serves to increase the 
scatter at high neutral fractions, when HII bubbles are small, as expected 
from the preceding argument. The inclusion of velocities is uncertain in any 
case because we ignore the possibility of galactic winds and correlations of 
the velocity field on large scales. 
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Figure 5. Scatter plot of the distance to the nearest neutral clump, Rb\ , as 
a function of the damping wing optical depth, for several different phases 
of reionization: Shi = 0.72, 0.51, 0.26 (top to bottom). Each panel was 
created using 1000 randomly chosen LOSs. The large scatter in Ri,i at a 
fixed to (z,s) [or in to (zg) at a fixed R^i] is one source of the difficulties 
with damping wing measurements. 

difficulty involved in accurately estimating xm in an inhomoge- 
neously ionized IGM. We have already seen that there is a deter- 
ministic and accurate mapping (Rbi, to) «-> xd = xm (c.f. eq.[4} 
in a uniformly ionized IGM. However, as discussed previously, un- 
der the more realistic picket-fence IGM absorption scenario, this 
mapping becomes stochastic. To illustrate this, in Figure [5] we 
show a scatter plot of the distance to the nearest neutral clump, 
Rbi, as a function of the damping wing optical depth. Each panel 
was created using 1000 randomly chosen LOSs in several different 
phases of reionization: xm = 0.72, 0.51, 0.26 (top to bottom). The 
large scatter in Rn at some fixed to(zs) is one source of difficul- 
ties with damping wing measurements: it implies that LOSs with 
identical apparent optical depths (at a specified wavelength) pass 
through very different patterns of HI. As the analytic models pre- 
dicted in the previous section, this scatter increases with decreasing 
xbi. 

5.1 The Distribution of the Inferred Neutral Fraction 

In Figure [6] we plot the probability distribution of 8 XD = 
(xo/xm — 1) for several different phases of reionization: xm = 
0.20, 0.26, 0.34, 0.42, 0.51, 0.61, 0.72 (right to left at large 8 XD ). 
The mean of 8 XD is greater than zero for all phases, as predicted by 
the toy models in the previous section: xn > xm- In addition, the 
bias and scatter increase as reionization progresses. Note also that 
the probability distributions of 8 XD are highly non-gaussian due to 
the restricted range < xd < 1. 

The small spike at 8 XD — 4 in the xhi = 0.2 distribution re- 
sults from LOSs containing a single neutral patch at the end of their 
path length, to which our prescription assigns xr> = 1. In reality, 
the optical depth along such LOSs is so small that to (and espe- 
cially Rbi) would probably not be measurable in the first place, and 
in any case our assumptions of fully ionized bubbles and a constant 
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Figure 6. Probability distribution of 5 XD = (xo/xm — 1) for several 
different phases of reionization: x m = 0.20, 0.26, 0.34, 0.42, 0.51, 0.61, 
0.72 (right to left at large S XD ). Note that the mean of 8 XD is non-zero, 
and that the distribution becomes wider and more biased as reionization 
progresses. 

xm along each geodesic break down in this regime. However such 
LOSs comprise less than 0.7% of the total at xm = 0.2, the small- 
est neutral fraction we study, and so do not significantly affect the 
statistical estimates below. We do note that the abundance of these 
LOSs does depen d on our semi-nu meric algorithm. In this example, 
the algorithm of IZahn et al J ( 120070 has 14% of LOSs in this regime 
at xm = 0.2. However, if these unusual LOSs (which are, as we 
have argued, probably useless for this measurement) are removed 
from the sample, both algorithms agree on the net bias. 

In Figure UJ we plot the bias expressed as (in — xm) (top 
panel), and the fractional scatter in xo (bottom panel). The solid 
curve is generated from all of the LOSs. This net bias is always 
positive and matches our toy model with exponentially distributed 
neutral patches fairly well. However, we do not see evidence for 
a turnover at small xm, and the simulation curve also increases 
somewhat more slowly then the toy model. 

The first of these differences has a simple explanation. Weak 
damping wing absorption might not be detectable with finite signal- 
to-noise observations (or it may not be separable from an uncertain 
source continuum), and the corresponding sources would likely be 
labeled as post-reionization objects. Thus, the dot-dashed curve is 
generated by setting xp = for LOSs with to(zs) < 0.1 (see 
Mesinger & Furlanetto 2007b for the total optical depth distribu- 
tions). Imposing a minimum value of to imposes a minimum on 
xd, so the bias starts decreasing at low xhi- As the neutral fraction 
decreases, the number of these LOSs increases rapidly; eventually, 
the apparent distribution will divide into a large set of apparently 
absorption-free spectra and a few spectra where the inferred neutral 
fraction is large. Of course, both sets must be taken into account to 
yield the strongest constraints. 

We also note that the absolute value of the bias is not particu- 
larly important, so long as it can be calibrated through simulations 
like this one. The crucial point will be to understand the sample 
and the model well enough to perform this calibration; otherwise 
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Figure 7. Top: The damping wing bias from assuming a constant neutral 
fraction throughout the IGM, expressed as (xjj— Shi), as in Fig.[T] Bottom: 
Fractional scatter in xo . In both panels, the solid curve is generated from all 
LOSs, the dashed (dotted) curve is generated from LOSs with to(zs) < 
5 (1), and the dot-dashed curve is generated assuming xo = for LOSs 
with tjj(zs) < 0.1. Note that the total mean bias is always positive, but 
removing parts of the distribution to mimic real observational data sets can 
drive the apparent mean bias to negative values at large Shi- 



systematic uncertainties will remain in the interpretation of the ob- 
servations. As noted above, the assumptions inherent in the particu- 
lar radiative transfer or semi-numerical algorithm used to generate 
the ionization field become increasingly important at xm < 0.2; 
thus we predict that damping wing measurements in a xui < 0.1 
universe will be very difficult to interpret. We have not extended 
our models to this regime because our semi-nu merical algorithm 
is no longer very robust (as the comparison to the Zah n et al]|2007l 
algorithm above shows) and because evolution across the line of 
sight will play an increasingly important role; a "light-cone" anal- 
ysis will probably be required for such a regime. 

On the other hand, real-world instruments are not infinitely 
sensitive, so samples generated by Lya emission line searches (e.g., 
narrow band surveys) might only contain objects whose optical 
depth is less than some maximum value following from the instru- 
mental sensitivity. To model this possibility, the dashed and dotted 
curves in Figure|7]are generated only from LOSs with to(zs) < 5 
and 1, respectively. Note that the apparent bias calculated from such 
truncated distributions becomes negative at large xbi, because im- 
posing a maximum to also imposes a maximum on xo (for a given 

Rbl)- 

The fractional scatter is comparable for all the curves. It in- 
creases most steeply at low values of xhi for the dot-dashed curve, 
because setting xo = for LOSs with weak absorption induces a 
bi-modal distribution of xo, with small values of xo shifting to a 
spike at =0 (c.f. Fig. [6]». Again, the scatter can be calibrated 
with these types of models, but it requires fairly large samples to 
interpret the red damping wing reliably. 
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Figure 8. Scatter plot of the power law index a from the damping wing 
profile parameterization in eq. jilt , fit using two points at z = zg and 
z = zs + 0.1. Panels assume Shi = 0.72, 0.51, 0.26 (top to bottom). 
The horizontal lines denote a = 1, as would be expected from a uniformly 
ionized IGM. 



5.2 The Damping Wing Profile 

In £]3.3r we noted that our toy models of picket-fence absorption pro- 
duce a steeper absorption profile than expected in a homogeneously 
ionized IGM. Here we confirm and further quantify this result us- 
ing our semi-numerical simulations. Specifically, we parametrize 
the absorption profile as (c.f. eq.|3]l: 

ro(z) oc R-f . (11) 

Note from eq. l[3} that in the homogenously ionized IGM, a = 1; 
however, from Fig. [3] we expect that during patchy reionization the 
mean power law index a is greater than unity. 

To test this with our simulations, we perform a simple two 
point power law fit to the profile shape, calculating to(zs) and 
to(zs + 0.1). A scatter plot of the resulting power law index, a, 
from 1000 randomly chosen LOSs is shown in Figure[8] Panels as- 
sume xm = 0.72, 0.51, 0.26 (top to bottom). It is obvious from the 
figure that indeed a > 1, with the profile steepening and the scat- 
ter increasing as reionization progresses. Note also that the mean 
and scatter of a change with xm, even at fixed to(zs)- Although 
LOSs can have similar integrated HI columns lengths at different 
epochs of reionization, the distribution of neutral hydrogen along 
this subset of LOSs does evolve. In general, LOSs intersect a fewer 
number of longer neutral patches at high Shi than LOSs with the 
same to at low xbi. Shorter neutral patches, especially those close 
to Rn, translate in turn into a steeper absorption profile (see the 
discussion in $3.3\ . 

Of course, the systematic variation of a with xm implies that 
the spectral shape can also be used to constrain the latter. However, 
because the dispersion in profiles is always at least as large as the 
differences in the means (even at fixed to), it will require large 
samples to take advantage of this information. 
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Figure 9. Probability distributions of 5 XD = (xo/xm — 1) generated 
from LOSs originating from halos with masses M = 2.6 X 10 11 , 2.5 X 10 10 , 
and 2.3 X 10 9 Mq (thick to thin curves), at several different phases of reion- 
ization: xm = 0.26, 0.51, 0.72 (right to left at large 8 XD ). 



5.3 Variation with Halo Mass 

iMesinger & Furlanettd d2007bl) showed that the mean and disper- 
sion of td(zs) were functions of halo mass at a fixed neutral frac- 
tion. Most of this variation is due to differences in the HII bubbles 
that these halos reside in (because larger halos are more clustered 
and so tend to sit in larger bubbles ; lFurlanetto et al.l2004) . But any 
excess difference would lead to another bias in the interpretation of 
damping wings. 

In Figure [9] we plot the probability distributions of 8 XD = 
(xd/xhi — 1) generated from LOSs originating in halos with 
masses M = 2.6x10", 2.5 xlO 10 , and 2.3 xlO 9 M (thick to 
thin curves), at several different phases of reionization: xhi = 0.26, 
0.51, 0.72 (right to left at large 8 XD ). The M = 2.6 x 10 11 M ha- 
los are the largest halos in our simulation at z = 9, with our 250 
Mpc box containing eight of them. 

There are clearly some differences in the inferred values of 
xd from the different types of halos: LOSs originating from more 
massive halos have somewhat narrower distributions, with smaller 
means, at fixed aim.. This is because more massive halos gener- 
ally sit inside larger bubbles (with larger Rbi) and so the Lya ab- 
sorption cross-section is flatter when photons enter their first neu- 
tral patch. As discussed previously, it is the varying cross-section 
that causes the bias and scatter in measurements; if the Lya cross- 
section were completely flat, xd would always equal xm- 

However, Figure [9] also shows that the distributions are much 
more sensitive to xm than to the halo mass and hence can be ro- 
bustly used to estimate xm from observational data sets even with 
little or no knowledge about the underlying halo. But the overall 
bias does have a non-negligible dependence on mass, so such in- 
formation will be useful. We quantify this in Figure [10] where the 
solid curves correspond to the same three mass scales as in Fig- 
ure^ M = 2.6 x 10 11 , 2.5 x 10 10 , and 2.3 x 10 9 M (thick to thin, 
or bottom to top). The mean bias, id- Shi, increases by a factor of 
~2 as the host mass scale is decreased from 2.6 x 10 to 2.3 x 10 9 



Figure 10. Damping wing bias statistics. Solid lines were created using all 
LOSs originating from halos with masses M = 2.6 X 10 11 , 2.5 X 10 10 , and 
2.3 xlO 9 Mq (thick to thin). Dashed and dotted lines were created using 
only LOSs originating from halos with masses M = 2.6X10 11 Mq and 
with Rf,i > 40 Mpc; the dotted lines additionally assume that xd = for 
LOSs with Tn(zs) < 0.01. Top: The damping wing bias from assuming 
a constant neutral fraction throughout the IGM, expressed in (xd — Shi)- 
Bottom: Fractional scatter in xd . 

Mq, and the scatter at small xm also decreases somewhat. Thus, 
if the properties of the host halo can be measured, it will certainly 
help to extract stronger constraints. In the following subsection, we 
examine such a special case. 

5.4 Quasars 

As m entioned in the introduction, IMesinger & Haimai] d2004l 
2007) already claim to have detected damping wings in two high- 
redshift quasars. Their model assumed a uniform UV background 
flux for the purposes of calculating the damping wing, so their final 
constraint is comparable in spirit to our xd parameter, though it is 
not clear if their results are dominated by the damping wing profile 
shape or the inferred td (and likewise xd), the later being partially 
degenerate with other free parameters in the analysis. 

Obviously, it would be interesting to study the effectiveness 
of such damping wing constraints when an inhomogeneously ion- 
ized IGM is included, even more so considering that they can be 
applied to future high-redshift data sets. Note that, although our 
semi-numerical simulation boxes are at zs = 9, well beyond the 
SDSS (and possibly future) quasars, the ionization topology and 
optical depth statistics are weak functions of redshif t in this range 
dMcOuimi et al.ll2007l;lMesinger & Furlanettoll2007bh . 

Unlike normal galaxies, which are the focus of most of this 
work, bright quasars lie in highly biased regions with correspond- 
ingly large Rti . We have already seen in Figures [9] and \10\ that 
this decreases the scatter and bias. In Figure [10] the thick curves 
show the bias and scatter computed from LOSs originating from the 
most massive halos in our simulation box, M = 2.6 X 10 11 Mq. 
For the present analysis, where the bright quasar necessarily pro- 
duces a large HII region, we are interested only in rare LOSs with 
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large Rti . To guarantee convergence of our measurements in these 
unusual cases, we extend our path length of integration to 400 
Mpc, although we find that this only has a noticeable effect for 
the xhi = 0.2 data point. The thick solid line was created using 
all LOSs. The dashed lines were created using the subset of LOSs 
with R bl > 40 Mpc, typical of the high-z SDSS quasars[f] 

Both the bias and the scatter decrease compared with the solid 
curve. Requiring that be large means that the Lya absorption 
cross-section at Rbi is flatter than usual, with (xd — Shi) smaller 
by ~ 0.05 throughout. We caution however that the curves in Fig- 
ure[l0]are calculated at zs', thus if one is estimating xd blueward 
of the line center (i.e. using td(z < zs)), the bias is likely to lie 
somewhere between the solid and dashed curves. Note also that the 
bias shown with the dashed curve becomes negative at xm > 0.6 
(see Fig.|7]and discussion thereof). The dotted lines in Fig.|10|were 
also created using LOSs with Rti > 40 Mpc, but with the addi- 
tional assumption that xd = for LOSs with td(zs) < 0.010 
Having an effective minimum td results in the same decrease of 
bias and increase in scatter as was seen in Figure [7] 



6 DISCUSSION 

In this paper, we have examined how the shape of the Lya red 
damping wing can be used to constrain the IGM before reionization 
is complete. In the past, it has usually been assumed that the absorb- 
ing gas can be well-approximated by a uniform density medium 
with constant ionized fraction. However, recent reionization mod- 
els have shown that ionized bubbles can be quite large, so the latter 
is not a good approximation. We have therefore critically examined 
how well the damping wing constrains the neutral fraction during 
inhomogeneous reionization. 

We have identified two major issues with its interpretation. 
First, there is substantial scatter in the optical depth along differ- 
ent lines of sight. Most of this is due to the scatter in the distance 
between the source and the nearest patch of neutral gas; however, 
there is still non-negligible scatter even if this distance can be mea- 
sured from the shape of the damping wing. In our semi-numeric 
simulations, the fractional r.m.s. fluctuation in xm thus estimated 
increases from 0.1 to 1 over the range 0.9 > xm > 0.2. Fortu- 
nately, this statistical uncertainty can be reduced simply by finding 
more lines of sight. 

The other problem is more severe: we have shown that the 
"picket-fence" absorption from inhomogeneous reionization adds 
a systematic, and often large, bias to measurements of the neutral 
fraction. Although the damping wing is indeed sensitive to a large 
path length through the IGM, it is most sensitive to the closest gas. 
As a result, simple fits to the damping wing will always overesti- 
mate the true neutral fraction in a partially ionized universe, with an 
error of ~ 30% near the middle of reionization. This bias is generic 
to any inhomogeneous model. The bias is reduced and can even be- 
come negative if observations only probe a subset of the entire halo 
population, such as quasars with large HII regions. 

Both the systematic and statistical uncertainty can be reduced 



5 LOSs with such large Rti are very rare at high xm and our box only 
contains them when xm < 0.75. 

e Note that we assume the damping wing is studied blueward of the Lya 
line center for quasars, so its footprint can be non-negligible even with a 
small line-center optical depth Td(zs). To(zg) ~ 0.01 would roughly be 
expected if xm ~ 0.1 and Rn ~ 40 Mpc, as in the observed systems. 



by a careful fit to the damping wing spectral profile, which is typ- 
ically steeper than the naively expected (AAobs) -1 profile. How- 
ever, because the absorption typically comes from many neutral 
patches, a large number of parameters are required for a detailed 
fit, and given the relatively modest difference from the expected 
behavior, these will be difficult to measure, probably only possi- 
ble in systems with intrinsically large optical depths. Moreover, the 
scatter in the profiles, even at fixed td, is sufficient that large sam- 
ples will be required to put strong constraints on reionization from 
the spectral shape. 

Of course, the bias and scatter also depend on xm and so can, 
at least in principle, be used to constrain it. For example, large dis- 
persion in the inferred neutral fractions could be an indicator of 
Shi <0.2. If an independent estimate of Shi exists, one could re- 
verse the direction of analysis, and use the bias and scatter to con- 
strain the reionization model and topology. 

Fortunately, for a given model of reionization, the dispersion 
and bias can be calibrated by theoretical models. We therefore ar- 
gue that the most efficient way to constrain reionization with the 
damping wing is through comparison with detailed models. Of 
course, any such constraints will be model-dependent, but we be- 
lieve that the morphol ogy of reionization is now sufficiently well - 
understood (see, e.g., Furl anetto et al.|[2006l : lMcOuinn et al ] l2007h 
that these uncertainties will likely not dominate the statistical un- 
certainties from the small number of accessible sources, at least 
in the relatively near future. For example, the reionization mor- 
phology is nearly ind ependent of redshift jFurlanetto et alj|2004l : 
McOu inn et alj 120071) . Also, we have found only a modest de- 
pendence of the xd distribution on halo mass (mostly due to the 
variation in bubble size with mass). However, toward the end of 
reionization, when the absorption is dominated by rare, narrow 
sheets of neutral hydrogen, the details of the radiative transfer al- 
gorithm (or an approximation to it, as in our models) and of the 
sample selection will be extremely important. Nevertheless, the 
task is challenging, as the damping wing profile must be separated 
from the rapidly varying resonan ce absorption for quasars (as in 
iMesinger & HaimarJ [2004 l2007h or from intrinsic absorbers for 
GRBs. Fortunately, in the latter case ~ 20% of moderate-redshift 
GRBs ha ve only modest ab sorbers and will still be useful for these 
purposes dChen et al"1l2007l) . 

So far, the damping wing analysis has been performed on 
three high-redshift quasars: Jl 148+525 1 (z s = 6.42), J1030+0524 
(zs = 6.28), Jl 623+31 12 (2 S = 6.22) ( IMesinger & HaimarJ 2004 . 
20071), as well as GRB 050904 (z s « 6.3) dKawai et aljboo" 



Totani etal 2006). This paper highlights the need to calibrate these 



and future damping wing analysis with simulations of the reioniza- 
tion morphology. Obviously we cannot set firm constraints without 
detailed simulations of the observations. Nevertheless, the mean 
bias we find from our simulations seems to work in the direction 
of stre ngthening the upper limit (on Shi) from the iTotani et"al] 
d2006h measurements, and weake ning the lower limit from the 
IMesinger & Haim an 1 2004. 120071) constraints at xm < 0.6 (al- 
though, interestingly, it would strengthen them if xhi > 0.6; see 
the sign change for the bias in Fig. |10t . Conversely, the steeper- 
than-expected absorption profile s eems to work in the direction of 
weakening the ITotani et al.1 d2006t) constraints (especially because 
it must be distin guished from strong i nterna l absor ption) while 
strengthening the Mesinger & Haiman] d2004l . 120071 ) constraints. 
The absorption profile might be more relevant than the bias for 
these studies, as an overall bias can be partially degenerate with 
other free parameters in the fit: that is, when the absorption profile 
can be detected to high precision, its shape will certainly be use- 
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ful in constraining xm- The scatter in both effects would probably 
somewhat erode the confidence contours for all of these studies. 
On the other hand, our model predicts large scatter between differ- 
ent LOSs at the end of reionization, which is consistent with the 
measurements at z ~ 6.3. More precise limits will require simulta- 
neous fits to the intrinsic absorption and the range of possible IGM 
absorber profiles, and we defer them to future work. 

Another intriguing possibility is to try to measure damping 
wing ch aracteristics from stack ed spectra of many Lya-emitting 
galaxies. McOu irm et alj d2007l) have shown that the wing shape 
is difficult to separate from uncertainties in the line for individ- 
ual objects, and the scatter we have described will also make the 
interpretation of individual faint emitters problematic. But, if the 
characteristics of the population are relatively constant, stacking 
may increase the signal to noise sufficiently to allow a detection of 
a "mean" damping wing at each redshift, even far redward of line 
center. 

SRF thanks Crystal Martin and Josh Bloom for conversations 
that stimulated this work. We thank Z. Haiman for helpful com- 
ments on this manuscript. This research was partially supported by 
grant NSF-AST-0607470. 
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